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ABSTRACT 

We report simulations of diffusive particle acceleration in oblique magnetohydrody- 
namical (MHD) shocks. These calculations are based on extension to oblique shocks of 
a numerical model for "thermal leakage" injection of particles at low energy into the 
cosmic-ray population. That technique, incorporated into a fully dynamical diffusion- 
convection formalism, was recently introduced for parallel shocks by Kang & Jones 
(1995). Here, we have compared results of time dependent numerical simulations using 
our technique with Monte Carlo simulations by Ellison, Baring & Jones 1995 and with in 
situ observations from the Ulysses spacecraft of oblique interplanetary shocks discussed 
by Baring et ai, (1995). Through the success of these comparisons we have demonstrated 
that our diffusion-convection method and injection techniques provide a practical tool 
to capture essential physics of the injection process and particle acceleration at oblique 
MHD shocks. 

In addition to the diffusion-convection simulations, we have included time depen- 
dent two-fluid simulations for a couple of the shocks to demonstrate the basic validity 
of that formalism in the oblique shock context. Using simple models for the two-fluid 
closure parameters based on test-particle considerations, we find good agreement with 
the dynamical properties of the more detailed diffusion- convection results. We empha- 
size, however, that such two-fluid results can be sensitive to the properties of these 
closure parameters when the flows are not truly steady. Furthermore, we emphasize 
through example how the validity of the two-fluid formalism does not necessarily mean 
that steady-state two-fluid models provide a reliable tool for predicting the efficiency of 
particle acceleration in real shocks. 

Subject headings: Cosmic-Rays — particle acceleration — magnetohydrodynamics 
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1. Introduction 

The theory of diffusive particle acceleration at col- 
lisionless shocks has been quite successful in explain- 
ing many aspects of the cosmic ray (CR) popula- 
tion. These include, for example, the nearly power- 
law spectrum of the CRs detected at the top of the 
atmosphere, the relation between the break in the 
power-law around the ~ 10^'* eV knee energy to the 
maximum energy of the CRs achievable in supernova 
remnants (SNRs), and, also the non-thermal, power- 
law electron populations deduced from the radio syn- 
chrotron observations of SNRs (Drury 1983; Bland- 
ford & Eichler 1987; Berezhko & Krymskii 1988). 
Although the primitive forms of the theory are very 
straightforward and robust, the microphysics is actu- 
ally complex, and there are potentially important sim- 
plifying assumptions built into the various versions 
of the theory. In practice nonlinear interactions be- 
tween the thermal plasma and the nonthermal (CR) 
plasma are also often likely to be essential. Further- 
more, there are certain key features, such as the pro- 
cesses that inject particles from the thermal plasma 
into the nonthermal plasma (hereafter, simply "injec- 
tion") that are not well understood. 

Over the past several years significant strides have 
been made in direct observational tests of diffusive 
acceleration theory and in comparisons between the- 
oretical models that are sometimes based on very 
different approaches. For parallel shocks, in which 
the ambient magnetic field is aligned with the shock 
normal, the applicability of diffusive shock theory 
is now fairly well established. Ellison and coUabo- 



ratprs have demonstrated 



between 



good agreemeni 
Mo'nte Carlo particle shock simulations and measure- 
ments at the earth's bow shock (Ellison, Mobius & 
Paschmann 1990), as well as between Monte Carlo 
and hybrid plasma shock simulation techniques (El- 
lison, et at, 1993). Recently Kang & Jones (1995, 
Paper I) demonstrated that "continuum transport" 
approaches based on the diffusion-convection equa- 
tion also provide good models for the same bow shock 
measurements and that in parallel shocks such con- 
tinuum models agree well with both types of parti- 
cle approaches. Paper I also demonstrated that the 
time-dependent two-fluid derivative of the diffusion- 
convection model works well as a dynamical model 
for these shocks, as long as the necessary closure pa- 
rameters are properly defined. 

Paper I, in addition, contained an important step 



for developing a physically-based injection model within 
the diffusive transport formalism. Particle methods 
have demonstrated injection to be an integral part 
of collisionless shock formation {e.g., Jones & Elli- 
son 1991 and references therein). Continuum mod- 
els for diffusive acceleration, on the other hand, have 
generally depended for practical reasons on an effec- 
tive separation of the particle population into distinct 
thermal and CR components with the former treated 
by fluid mechanical techniques and the latter by ap- 
proximate plasma kinetic equation techniques or the 
energy moment of that equation {i.e., the two-fluid 
model). Paper I introduced a hybrid model of the 
standard continuum approach. It creates a virtual 
"injection pool" of particles that are neither fully 
thermal nor fully "CR", but represent the particles 
leaking out of the thermal population into the CR 
component. We demonstrated there that this contin- 
uum "thermal leakage" injection model could produce 
good agreement with simulations done using particle 
techniques and with direct bow shock measurements. 
We focus our discussion, by the way, on the ionic 
particle population, rather than electrons, since ions 
carry most of the momentum flux and seem to cap- 
ture the greater share of energy in the shock transi- 
tion, especially among the suprathermal, high energy 
population. To simplify discussion, we deal only with 
protons, although the methods being used can also be 
applied to other ions. 

The situation regarding oblique shocks is not so 
well studied and technically more difflcult to model. 
Quasi-parallel and quasi-perpendicular collisionless 
shock structures may be quite different {e.g., KennelJ 
Edmiston & Hada 1985). There are complications 



in the details of diffusive particle propagation (such 
as anisotropic diffusion with respect to the magnetic 
fleld) and acceleratio n that derive from the magnetic 
fleld geometry {e.g., Jokipii 1987 ). In addition, the 
magnetic fleld can apparently modify and reduce in 
nonlinear ways the amount of energy that CRs are 
able to extract from the flow through highly oblique 
shocks {e.g., Webb, Drury & Volk 1986; Jun, Clarke 
& Norman 1994; Frank, Jones & Ryu 1994, 1995). 
As a signiflcant step towards a better understand- 
ing of the acceleration physics in oblique shocks Bar- 
ing, et al, 1995 (BOEF95 hereafter) have extended 
the Monte Carlo studies to oblique shocks through 
test-particle simulations in which a gyro-orbit compu- 
tation was adopted and large-angle scatterings were 
assumed for the particles. They demonstrated that 
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Monte Carlo shock models with reasonable scattering 
properties can match in situ observations of oblique 
interplanetary shocks from the Ulysses spacecraft. El- 
lison, Baring, & Jones 1995 (EBJ95 hereafter) have 
applied the same techniques to study the accelera- 
tion rates and injection efficiencies in oblique shocks. 
They confirmed earlier findings from different meth- 
ods that, for anisotropic diffusion, the acceleration 
rate for individual particles increases with the mag- 
netic field obliquity {e.g., Jokipii 1987; Frank, Jones 
& Ryu 1995). Thus, quasi-perpendicular shocks may 
be capable of generating higher energy particles than 
quasi-parallel shocks in a given time, even though 
their net energy extraction efficiencies may be re- 
duced by a strong field. In their Monte-Carlo simula- 
tions EBJ95 also saw that the efficiency of "thermal- 
leakage" injection decreases with the obliquity, mak- 
ing it harder for quasi-perpendicular shocks to gen- 
erate seed CR particles. They found that the injec- 
tion rate depends upon the Mach number, field obliq- 
uity angle, and strength of Alfven turbulence respon- 
sible for scattering. The last of these enters, because 
it controls the amount of cross-field diffusion, which 
in quasi-perpendicular shocks becomes necessary for 
particles to escape into the upstream region in order 
to be accelerated. A recent complementary discussion 
of the characteristics of electron injection at quasi- 
perpendicular shocks has been provided by Levinson 
(1996). 

The present discussion extends the analysis of Pa- 
per I to oblique shocks. We apply our thermal leak- 
age injection model to oblique shocks and make a 
preliminary comparison with the behaviors reported 
by EBJ95. In addition we turn the model to the 
same interplanetary shock measurements presented 
by BOEF95. These results will provide a useful foun- 
dation for future studies of oblique shock physics 
based on continuum models, and expand on a pre- 
liminary report made earlier (Jones & Kang 1995). 
The plan of the paper is this. In §2 we summarize the 
comparison between particle and continuum models 
and how that is important to understanding injec- 
tion. Section 3 outlines our numerical methods, while 
§4 presents our results, and §5 provides a summary 
and conclusion. 

2. Continuum Models and Particle Injection 

In continuum treatments of diffusive shock ac- 
celeration, the diffusion-convection equation (equa- 



tion 3-1; to be discussed in §3) for the particle mo- 
mentum distribution is solved along with the dy- 
namic equations for the underlying plasma. They 
have some distinct practical advantages over Monte 
Carlo methods, especially in time dependent prob- 
lems and those that involve complex, or multidimen- 
sional flows. There are well developed, robust and 
relatively inexpensive continuum computational tech- 
niques available that can be applied very flexibly, 
for example. There is a growing literature based on 
time dependent diffusion-convection equation treat- 
ments of quasi-parallel shocks {e.g., Falle & Giddings 
1987; Kang & Jones 1991; Kang, Jones & Ryu 1992; 
Duffy, Drury & Volk 1994). The simpler and much 
more economical, two-fluid derivative of the diffusion- 
convection method has seen even greater application 
in time dependent dynamical problems, because it is 
practical to use in many complicated situations where 
suitable numerical gasdynamics methods are suitable 
{e.g., Drury & Falle 1986; Dorfi 1990, 1991; Jones 
& Kang 1990, 1992, 1993; Ryu, Kang & Jones 1993; 
Jones 1993; Jones, Kang & Tregilfis 1994). Recently, 
both of these continuum techniques have been ex- 
tended to time dependent MHD and oblique shocks 
(Frank, Jones & Ryu 1994, 1995). Jun, Clarke & Nor- 
man (1994) also reported two-fluid results for perpen- 
dicular MHD shocks. Because of the wide-spread use 
of two-fluid methods their validity and limits of appli- 
cability take on a particular importance. Paper I as 
well as some other earlier papers {e.g., Duffy, Drury 
& Volk 1994) have demonstrated basic validity of the 
model in quasi-parallel shocks. There are some im- 
portant caveats, however, as we shall discuss later. 

The diffusion-convection equation is based on the 
assumption that a particle momentum distribution 
is kept almost isotropic in the local fluid frame by 
scattering. When particle velocities are much greater 
than the bulk speed of the plasma motion and scat- 
tering is efficient these requirements are very reason- 
able. It is reassuring, but perhaps not terribly sur- 
prising that statistical particle approaches like Monte 
Carlo and continuum (diffusion-convection) methods 
agree with each other and with real data in that case. 
In the opposite limit of random particle speeds less 
than those associated with bulk motion relative to 
the shock, the theoretical situation is again tractable. 
Since charged particles generally are scattered more 
rapidly at small velocities, we can expect them to 
be "thermalized" effectively and isotropized with re- 
spect to a well-defined mean motion. Then continuum 
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fluid dynamical computational methods should work 
well, and statistical methods should be convergent 
with them. The situation is much more problematic 
at "intermediate" velocities where particle streaming 
may be large, and mean free paths are too long to 
allow scattering to relax the distribution to a thermal 
form. (Particles sample a range of environments, so 
that no simple thermal equilibrium is appropriate.) 
This is the arena of injection, since it includes parti- 
cles that are capable of becoming CRs through multi- 
ple shock crossings. The reality of injection in shocks 
is not much in doubt {e.g., Jones & Ellison 1991), 
but despite recent theoretical strides {e.g., Malkov & 
Volk 1995), the details of injection remain beyond 
straightforward models. This problem is difficult be- 
cause nonlinear interactions between particles, reso- 
nant hydromagnetic waves and the underlying plasma 
associated with the shock formation process itself are 
very complex and yet to be deciphered fully. In that 
context we seek now only a simple, functional model, 
but one that captures essential physics of the process. 
Previous injection schemes within the continuum for- 
malism have generally been based on ad hoc assump- 
tions that a fixed fraction of the kinetic energy flux 
or the total particle flux through the shock are trans- 
ferred from the "thermal" to the "nonthermal" pop- 
ulations, so our new approach represents a clear step 
towards reality. There is considerable value in devel- 
oping a serviceable, but physically based model for 
injection within the continuum transport paradigm. 

The thermal leakage injection model introduced in 
Paper I is conceptually simple and represents only 
a small change from previous diffusion-convection 
methods. As before we simultaneously solve the cou- 
pled diffusion-convection/MHD equations. However, 
in this new technique we follow the entire proton 
momentum distribution with the diffusion-convection 
equation, but continuously redistribute the particles 
at low momenta into a thermal distribution according 
to the pressure and density solution from the MHD 
equations. That introduces a population of diffusive 
particles at intermediate momenta between the ther- 
mal particles and those properly termed CR particles. 
Since they are diffusive those particles can sample the 
fluid velocity on both sides of the shock, if they are 
given scattering properties suitably matched to the 
numerical shock thickness. Those intermediate mo- 
mentum particles gain energy in a manner that re- 
sembles what happens to CRs, but their distribution 
is directly matched onto the thermal distribution, as 



it physically must be. The injection efliciency is de- 
termined by the momentum at which this "injection 
pool" distribution is matched to the thermal distri- 
bution. Eventually, we hope to be able to provide an 
independent model for this matching, but for now, it 
can be sufficient to show that the model successfully 
reproduces important physical behaviors for physi- 
cally reasonable matching conditions. 

3. Model Description 

3.1. Numerical Methods 

We follow evolution of the particle distribution 
with the standard diffusion-convection equation {e.g., Parker 
1965; 



Skilhng 1975| ; [Jokipu 1982| ), 

V-(kV/)-u„ 



. Of 



•V/, (3-1) 



where f{x,p,t) is the isotropic part of the distri- 
bution function measured in the converted frame, 
u + Uw, while d/dt is the total time derivative in 
the ffuid frame, u. The propagation of scattering 
centers relative to the plasma is represented by u.u]. 
Generally, the scattering centers are assumed to be 
Alfven waves resonant with the particles, so repre- 
sents the center-of- momentum motion of those Alfven 
waves. Our simulations assume planar symmetry, 
so K = K,{p)i termed the diffusion coefficient, is the 
projection of the spatial diffusion tensor onto the 
shock normal. That direction is taken to be paral- 
lel to the X axis. Throughout the paper we express 
momentum, p, in units of the proton rest mass en- 
ergy, m,(? jc {— 9.38 X lO^keV/c), and the distribu- 
tion function / in units of the particle number den- 
sity, so that 47r J / dp — p/m. We solve equation 
3-1 on an Eulerian grid using a second-order com- 
bined Lagrangian Crank-Nicholson/monotone-remap 
scheme whose details can be found in Kang & Jones 
(1991) and Frank et al, (1995), respectively. The 
distribution f{p) is supposed to exist over a range 
Po l£ P l£ P3 that includes both the thermal distri- 
bution and the CR distribution. The thermal distri- 
bution expressed in terms of g{p) = p'^.f{p), which 
measures the partial pressure, dP/dp, has its max- 

/ir 



imum at pth = V 4r = y^Am ks T/mc, where the 
gas temperature T is expressed in units mc^ /ks. We 
will identify those particles dynamically as CRs that 
satisfy p> P2 >> Pth, and will establish p2 below. 

The dynamics of the underlying plasma is fol- 
lowed by an explicit, second-order accurate MHD 
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code based on a conservative up-winded, Total Varia- 
tion Diminishing (TVD) scheme (Ryu & Jones 1995) 
that has been modified to include the dynamic ef- 
fects of the CR pressure (Frank et ai, 1995). Read- 
ers are referred to their papers for the basic MHD 
equations and the detailed description of the numer- 
ical method. Based on a linear Riemann solver used 
to compute "up-winded" mass, momentum and en- 
ergy fluxes at zone boundaries, the code generally 
captures cleanly all the families of MHD discontinu- 
ities. Strong shocks are usually contained within 2 
to 3 zones, and other discontinuities within a slightly 
broader space (~ 4 — 10 zones, depending on the fea- 
ture). The code is conservative in the sense that it 
maintains exact net fluxes through the grid to ma- 
chine accuracy. The TVD label refers to the manner 
in which the code avoids introducing physically spu- 
rious oscillations by preserving monotonicity in phys- 
ical flow variables through discontinuities. 

There is one addition to the code discussed in 
Frank et ai, (1995); namely, "Alfven Wave Trans- 
port" (AWT) terms, as represented in equation 3-1 
by Uyj. Those are handled in the same way as dis- 
cussed by Jones (1993) and Paper I for parallel shocks, 
with the proviso that H^,, aligns with the local mag- 
netic field vector. Additional AWT terms provide for 
gas heating due to dissipation of the energy trans- 
ferred from CRs (the last term in equation 3-1) to 
Alfven waves (thence to the plasma) and also trans- 
port of the energy and momentum content within the 
waves. In the present simulations we have neglected 
the energy and momentum carried explicitly within 
the wave field, but have included the energy and mo- 
mentum passed through the wave field (see Jones 1993 
for details). The magnetic field in this problem lies 
within a single plane containing the shock normal di- 
rection, X. So, without loss of generality we can define 
the magnetic field to be within the x — z plane. 

For oblique MHD shocks the diffusion coefficient 
takes the standard form 



'II' 



(3-2) 



where || and _L refer to diffusion along and across 
the magnetic field direction, respectively, and 6 = 
a.rcta.n{Bz/ Bx). Following Jokipii (1987), we assume 
a parallel diffusion coefficient of the form K|| = ^A||W, 
with the scattering length, A|| = Nrg, where is the 
gyro-radius of a particle and v is its speed. Then, 
from standard kinetic theory, the ratio of the paral- 
lel to perpendicular components is determined by the 



ratio iV > 1 as 



(1 



(3-3) 



Equation 3-2 can be rewritten as 



K=[N cos^ e + ( ) sin^ e] KB , (3-4) 
1 + iV-^ 



where kb — 5^9"^ is the Bohm diffusion coefficient. 
The limit ^ corresponds to Bohm diffusion, 
where kj_ /^n- Cross-field diffusion is determined 
in this model by the strength of Alfvenic turbulence, 
since, N ^ Eb / {kE^k), where Eb and E^k are the 
total energy density in magnetic fields and the Alfven 
wave energy density at the resonant wave number, 
fc, respectively. When scattering is weak, so that 
N » \, there is little cross-field diffusion, whereas 
strong scattering leads to cross field diffusion compa- 
rable to field-aligned diffusion. If is a constant, it 
follows for nonrelativistic particles that k oc . 

3.2. A Numerical Injection Model 

As we stressed above, the detailed physics of the 
injection process is not yet well understood, and dif- 
fusive transport models cannot, by themselves, ac- 
curately treat the particles directly involved in the 
process. So, as a practical approximation we assume 
a simple but reasonable scenario in which a small 
population of near-thermal particles gain excess en- 
ergy via interactions with resonant waves and form 
a suprathermal tail on the Maxwellian distribution 
in the vicinity of the shock front. They provide the 
seed particles injected into the CR population. As 
noted before, in our model the particle distribution 
over the full range of momenta including the thermal 
plasma is followed explicitly. Below a certain mo- 
mentum (in units of mc), pi = c\ pth — ci \/(4T) = 
ci\/AmkBT/mc, chosen high enough to include most 
of the postshock thermal population, the distribution 
is forced to maintain a Maxwellian form consistent 
with the local gas density and pressure determined 
from the MHD equations. Above pi particles are al- 
lowed to evolve according to the diffusion-convection 
equation, while only for p > P2 > Pi are they con- 
sidered dynamically as CRs. Particles between pi 
andp2, thus, constitute "candidate" CR particles, be- 
cause they are not locally thermalized. They can be 
injected into the CR population by crossing a mo- 
mentum boundary ai p — p2 through flow compres- 
sion. Below pi particles are compressed adiabatically 
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(rcvcrsibly) by the flow, except within the shock dis- 
continuity, where the shock jump conditions demand 
that the compression be irreversible. Above pi, on 
the other hand, compression leads to irreversible en- 
ergy changes in the particles, because diffusion is ir- 
reversible. This combination of effects is the source 
of energy for the particle acceleration, of course. 

We emphasize that the distribution function for 
thermal particles is used only to provide the refer- 
ence population needed to match onto the interme- 
diate population, while the thermal pressure, Pg, is 
included and handled through the MHD equations. 
On the other hand, the intermediate population only 
provides the seed particles for CR particles and has 
no dynamical effects on the flow in our method. Al- 
though the particle distribution is continuous over the 
full range of momenta, in continuum treatments short 
of a full solution to the Boltzmann equation one needs 
to separate the pressure due to thermal particles (Pg) 
and that due to CR particles (Pc), because their dy- 
namical behaviors are different. This, of course, ne- 
cessitates a definition for the CR population. We have 
done that by choosing p2 as the arbitrary boundary. 
At early stages of acceleration, when the particle dis- 
tribution is almost Maxwellian, the small CR pres- 
sure is sensitively dependent upon the chosen value 
of p2, but Pc is dynamically insignificant then. On 
the other hand. When Pg becomcis large enough to 
be important, it becomes much less dependent upon 
P2, because particles of much higher energy dominate 
the CR pressure. As a result, our calculations are not 
critically dependent upon the parameter p2. So, it 
is most convenient numerically to fix the value of P2 
--^ (3 — ^)pth,i, where pth,i is the thermal peak mo- 
mentum of the postshock gas of the initial pure gas- 
dynamic shock. The particle supply in the interme- 
diate momentum pool is sensitively controlled by the 
parameter pi, since f{pi) is part of the exponential 
tail of the postshock Maxwellian distribution. Ac- 
cording to comparison tests with measurements of a 
parallel occurrence of the earth's bow shock and with 
shocks computed by "particle" methods (see Paper 
I), appropriate values of the related scaling param- 
eter, ci = pi/pth, fall in the very reasonable range 
ci = 1.5 — 2. Thus the value of pi varies in time and 
space along with the local gas temperature. 

The model further requires us to match the numer- 
ical shock thickness, 6x, to particle scattering proper- 
ties, since the numerically realized injection rate will 
depend upon the ratio \{p\)/5x (see Paper I). Above 



pi particles are formally diffusive, but unless the scat- 
tering lengths of these particles projected onto the 
shock normal, A(p)cos^, exceed 5x, they cannot be 
effectively accelerated by the Fermi process. On the 
other hand thermal particles should not be able to 
cross the shock within a projected scattering length, 
since they should then not form into a Maxwellian dis- 
tribution. Hence, the numerical shock must be thicker 
than the projected scattering length of thermal parti- 
cles, but thinner than the projected scattering length 
of CRs. The structure and thickness of real shocks will 
be dependent upon the details of the strength and ge- 
ometry of the field, degree of turbulence, the strength 
of the shock, for example. That issue is beyond the 
scope of this study. The specifics adopted for the nu- 
merical shock thickness will be given later for each 
case. This "thermal leakage" type injection model is 
rather simple, but, according to the results reported 
in Paper I, apparently able to capture essential char- 
acteristics of real injection processes, provided that 
one makes a reasonable choice for the free parameter 
Pi- 

3.3. Initial and Boundary Conditions 

For our simulations the initial flow is specifled by 

a simple discontinuous MHD shock using standard 
jump conditions, which can be found from MHD 
Riemann solutions, for example, (see Ryu & Jones 
1995). The shock faces to the right, so that veloc- 
ities along the x axis are negative when the shock 
is nearly at rest in the grid. All of our shocks 
start at rest in the grid, but those developing dy- 
namically significant CR pressure become temporar- 
ily "over-compressed", as expected, causing them to 
drift slowly to the left. Three fluid parameters are 
needed to define the shocks. Those can be the sonic 
Mach number. Mi = Uix/csi, the strength of the up- 
stream magnetic field, Bi and the upstream obliq- 
uity of the magnetic field, 9i = arctanSia^/Bi^. The 
sound speed is c^i = ■^'yPg/p, where Pg is the gas 
pressure, p is the gas density and 7 is the gas adiabatic 
index, taken to be |. We also assume that initially 
the particle distribution function, f{p), is Maxwellian 
everywhere, with a temperature, T = {Pgin) / {pkB)- 
Thus, there are no CRs initially. They are injected 
through thermal leakage as part of the process of 
shock evolution. The above definition of the temper- 
ature, which was used in both EBJ95 and BOEF95 
implies that the electron pressure is negligible com- 
pared to the proton pressure. More recently. Baring et 
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aL, 1996 (BOEF96) have recomputed the properties 
of the Ulysses-observed shocks, including finite elec- 
tron pressure. Although that has changed some of the 
shock parameters, it should not alter any of our con- 
clusions. The influence of a finite electron pressure is 
certainly straightforward to include when warranted. 

The MHD variables are assumed to be continu- 
ous across the left and right boundaries of the spatial 
grid. This is a good assumption, since the shock is ap- 
proximately at rest in the middle of the grid, keeping 
any gradients in flow variables small near the bound- 
aries. The particle distribution, f{x,p), is also as- 
sumed to be continuous across the boundaries, which 
means diffusive particle fluxes vanish there. This no- 
flux boundary condition is numerically simple and ro- 
bust for the diffusion-convection equation. It remains 
reasonable as long as the particles are confined near 
the shock and away from the boundaries. 

The boundary condition for f{p) just below po is 
not relevant here, since the distribution is redefined 
continuously by the Maxwellian function at each time 
step for p < pi. At the highest momentum boundary, 
we assume f{p) = fips) for p > p:^. This condition is 
not very crucial either, since the divergence of the flow 
is rare around the shocks in plane-parallel geometry 
considered here. 

4. Results 

4.1. Comparison with Monte Carlo Simula- 
tions 

EBJ95 have calculated, by test-particle Monte Carlo 
simulations, the efflciency of injection at oblique shocks 
as a function of Mach number. Mi; field obliquity, Oi; 
and the degree of cross-field diffusion (as measured 
by N = X/rg). They found the injection to be more 
efficient for lower Mach numbers, for smaller obliqui- 
ties and for stronger cross-field diffusion (i.e., smaller 
N). In their Fig. 5 they showed the downstream in- 
tegral density distribution for particles accelerated in 
strong shocks for a range of obliquity. This informa- 
tion can be compared directly with the particle distri- 
bution functions of our simulations. Thus, we chose 
these shocks as the comparison models and found the 
value of ci for each value of 9i that gives the best 
fit to their results. The common shock parameters 
are Ui^ = 500 km s-\ Mi = 100, N = X/rg = 100, 
and Bi = 10^^ Gauss. This represents a very strong 
shock in the limit of weak cross-field diffusion and 
weak magnetic field. The obliquity values considered 



are 6*1 = 0°, 20°, 30°, and 35°; so that all arc "quasi- 
parallel" shocks. Larger obliquities were not consid- 
ered by EBJ95 for this shock system, since CR injec- 
tion was fomid to be completely suppressed. 

The EBJ95 simulations were test-particle, so for 
this comparison test only, we turned off the dynamic 
evolution of the fiow and kept the shock structure 
as the initial discontinuous jump (thus, with no CR 
pressure feedback, even though our code is designed 
to include fully the dynamical contributions of the 
CR). For these test-particle runs the shock thickness 
is effectively one grid cell. Since the shock thick- 
ness should be of order the mean scattering length of 
the postshock thermal momentum, pth, we adjust the 
grid spacing to be this length (e.g., Ax = X{pth) = 
Nrg{pth))- AH physical lengths in this problem scale 
with X{pth), so this model for the shock thickness will 
make the injection process scale with N. Our compar- 
isons with EBJ95 were with Monte Carlo simulations 
that did not include AWT, so we turned those effects 
off for this particular set of continuum transport sim- 
ulations. 

Monte Carlo simulations intrinsically consider a 
steady state, while our calculations are time-dependent 
Fig. 5 in EBJ95 shows the Monte Carlo, integral den- 
sity distribution up to 1000 keV. In order to make 
a good comparison we should integrate our simula- 
tions for a time comparable to that required to accel- 
erate a thermal particle to E > 1000 keV. In prac- 
tice, however, these become fairly expensive for cases 
with smaller obliquity, 0\, because the integration 
time is longer [tacc oc k), and so a greater spatial 
length is required to keep the CR particle distribu- 
tion small at the boundaries. Thus we evolved each 
shock for a time needed to accelerate particles to E ^ 
a few X 100 keV. That corresponds to (t/lO^s) = 
12, 8, 6, and 4 for 91=0°, 20°, 30°, and 35°, respec- 
tively. Fig. 1 shows the resulting particle distribu- 
tions at the shock position for the times and values 
of 6i specified above. The cases shown here are the 
models with ci chosen to match the results of EBJ95. 
The top panel shows p*f{xs,p), where Xs is the shock 
position. The canonical test-particle distribution is a 
power-law; in this strong shock case, f{p) oc p^. The 
distributions shown in Fig. 1, however, are some- 
what steeper than this canonical power-law, because 
the computed momentum range is finite and because 
the computed time interval does not a real steady 
state to be achieved. The slopes of pow-law fits of 
these distributions at p ~ 1.5 x 10~^, for example, are 
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q - 4.06 - 4.1, where f{p) oc jr^ . Following EBJ95, 
the bottom panel shows the integral of the distribu- 
tion function above a given kinetic energy, expressed 
in units of keV; namely, n(> E) = An p''^f{p')dp' , 

where E = mc'^{^p'^ + 1 — 1). The bottom panel also 
includes the analogous results reported by EBJ95 for 
their simulations. 

As mentioned before, our spectra start to cut off 
above E ^ a few x 100 keV, due to limited evo- 
lution time, while the steady-state, EBJ95 Monte 
Carlo results extend to higher kinetic energy values. 
More recently Ellison, Baring & Jones 1996 (EBJ96) 
have extended their test-particle simulations to fully 
dynamical Monte Carlo simulations. In those sim- 
ulations they included a "Free Escape Boundary" 
(FEB), which removes particles that propagate "too 
far" upstream from the shock. That preferentially re- 
moves the highest energy particles, since they have 
the longest scattering lengths. The net result is an 
energy ciitoff that qualitatively resembles the finite- 
time cutoff observed in the distributions we show in 
Fig. 1. Notice that the high energy side of the EBJ95 
"quasi-thermal" distributions cut off more sharply 
than Maxwellian. This presumably results from the 
rapidly increasing rate of thermal particle "leakage" 
with momentum in the Monte Carlo simulations. Our 
distribution, on the other hand, is not allowed to devi- 
ate from the Maxwellian form below p\ corresponding 
to £" ^ 2 keV, and we simply match the nonthermal 
distribution to it. But we see that within an energy 
factor of 2 or 3 of the thermal peak (pth ~ 1-5 x 10~^, 
Eth ~ 1 keV) our distribution converges fairly well 
to that found by EBJ95, below the cutoff imposed 
by finite acceleration time. Thus, on the whole our 
model shows itself to be a reasonable way to mimic 
the injection and acceleration processes. It produces 
a consistent particle spectrum at energies higher than 
thermal energies, in agreement with the Monte Carlo 
simulations, even though the details of the injection 
of suprathermal particles are not included. 

The values of ci adopted for the test-particle simu- 
lations that fit best with EBJ95 results are 1.4, 1.65, 
2.0, and 2.3 for i9i = 0°, 20°, 30°, and 35°, respec- 
tively. The increasing values of ci for higher obliq- 
uity are required to reduce the injection rates for those 
shocks. From Fig. 1 or from EB J95 Figs 5 & 6 it is ap- 
parent that the injection rate decreases by about two 
orders of magnitude between 9i = 0° and 9i = 35° for 
this Mach number and N value. In fact, EBJ95 argue 
within the test-particle picture that above 6i 30°, 



injection within strong shocks may be completely sup- 
pressed in the absence of cross-field diffusion. The 
reason for the obliquity dependence is that particles 
propagate along field lines until they scatter, except 
for a drift along the shock plane that can be elimi- 
nated by referring to the so-called de Hoffmann- Teller 
frame. In strong, oblique shocks tan 02 ~ r tan 9i , 
where 62 is the downstream field obliquity and r is 
the shock compression ratio; that is, the downstream 
obliquity is greater than the upstream obliquity and 
downstream particle motions are more nearly along 
the shock plane. Thus, as the obliquity increases, a 
relatively larger total particle speed after an initial 
scattering is needed to enable a particle to "swim up- 
stream" fast enough to re-cross the shock from down- 
stream. This tendency reduces the number of parti- 
cles available for injection (Baring, Ellison & Jones 
1994). In our simulations the same effect is estab- 
lished by higher values of ci for higher obliquity. As 
ci increases thermal leakage is reduced, because the 
number of particles in the injection pool is reduced. 
From the above explanation it is clear in this model 
that injection is less sensitive to obliquity when the 
Mach number is smaller or when the scattering is 
stronger (that is N is smaller) . EB J95 found in those 
situations that the injection rate is also greater. This 
implies that smaller values of ci should be chosen in 
our model for smaller M and for smaller TV, since the 
injection efficiency is mostly controlled by the value of 
Ci. However, we do not attempt here to find a quan- 
titative dependence of ci on M or A^, since the in- 
formation presented in EBJ95 is insufficient for that. 
Also, the best-fit values of Ci could vary with the nu- 
merical shock thickness. We leave for the future a 
more detailed analysis of these model properties. 

The above simulations, both ours and those in 
EBJ95, were of a test-particle character. On the other 
hand, it is clear that the energy represented in the 
super-thermal particle distributions is a substantial 
fraction of the total. Thus, test-particle results are 
not very meaningful as a measure of the properties of 
real shocks of this kind. This is not surprising, since 
previous studies of strong gasdynamic CR shocks have 
found them to be very efficient at transferring energy 
from the flow to CRs {e.g., Drury & Volk 1981). 

To gain some insights into the properties of these 
shocks when they are constructed self-consistently, 
we repeated our simulations, but with the fully dy- 
namic version of our MHD/diffusion-convection code. 
To keep the tests simple we used the same val- 
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lies of ci to model injection as in the test-particle 
simulations. For an obliquity less than 20°, how- 
ever, more than 10% of the particle number den- 
sity is within the CR population (see Fig.l) by the 
end of the simulation interval, using the injection 
rates found by EBJ95. So, one more assumption 
of the MHD/diffusion-convcction approach is invalid; 
namely, that the inertia within the high energy, CR 
population can be neglected. For this reason we have 
done dynamic test runs only for 9i =30°. For conve- 
nience in later discussions we refer to this model shock 
as EBJ95-D. Unlike the test-particle runs where the 
shock is one cell thick, strong shocks in fully dynamic 
runs are captured within about two cells in our code. 
Thus, for these tests the grid is adjusted so that a 
cell has thickness. Ax = 0.5X{pth) in order again to 
match the numerical shock thickness to the scattering 
length for thermal particles; i.e., 5x w X{pth)- 

Fig. 2 shows the flow structure around the EBJ95- 
D shock at tl{l(fs) = 4, 8, and 12, along with the 
initial MHD shock jump at t = 0. The frame of ref- 
erence is chosen so that the shock is at rest with- 
out CR modification to the flow. The physical vari- 
ables are expressed for simplicity in units of the fol- 
lowing normalization constants: Lq = 5 x 10^^ cm, 
po = 1.67 X 10~2'*g cm-3, Mo = 5 X lO^km s~\ and 
Pgo = 4.175 X 10~^erg cm~^. The numerical grid ex- 
tends from 0.0 to 3.0 in units of Lq. Only the region 
between 0.0 and 2.0 is shown in the figure, however. 
The assumed value of ci = 2.0, is the same as for the 
test-particle run. Predictably, the CR pressure is dy- 
namically important for this shock, so that it has a 
clear precursor. Similarly, the maximum compression 
is 4.6 instead of the test-particle value of 4 and the 
postshock gas pressure is lower than that of the initial 
shock. Note at i/(10^s) = 12, that the shock struc- 
ture is still evolving rapidly and the postshock CR 
pressure is already about 30 % of the gas pressure. 
Thus, it is clear that the test particle approximation 
is not valid even for this obliquity. Since the shock 
structure has been modified by the CR pressure from 
the initial shock jump, the shock is moving slowly to 
the left from the initial position. As time goes on, the 
modified structure extends downstream (to the left in 
Fig. 2) due to advection, while the precursor in the 
velocity and the CR pressure extends upstream via 
the diffusion of highest energy particles. 

Fig. 3 provides a comparison of the particle distri- 
bution and the integrated particle density for EB.J95- 
D, as well as the test-particle calculation shown in Fig. 



2 with the same initial conditions. We first note that 
the gas temperature is lower, so the peak momentum 
of the Maxwellian distribution, pt^, is lower in the 
dynamic run, as we expect from the lower postshock 
gas pressure and higher compression shown in Fig. 2. 
The postshock gas is colder in the dynamic run and so 
the particles in the thermal tail have smaller r^, while 
the shock numerical thickness is the same length in 
both runs. Thus they are less likely to be able to cross 
the shock in the dynamic run than in the test particle 
run. This will reduce the injection rate in the fully 
dynamic run. Therefore, we expect for similar reasons 
that the injection rates in fully dynamic Monte Carlo 
simulations would decrease from those given in the 
EBJ95 test particle simulations, especially for small 
obliquities. 

4.2. Comparison with Ulysses Observations 

BOEF95 have compared proton distributions mea- 
sured directly by the Ulysses spacecraft at oblique 
interplanetary shocks with results from Monte Carlo 
simulations of similar shocks. The BOEF95 simula- 
tions are also test-particle ones. We have adopted 
the same shock parameters as they obtained, and cal- 
culated the time-dependent evolution of the particle 
distribution functions. Our runs include fully the 
dynamical feedback of CRs on the shock structure; 
however, because we are comparing our results with 
the Ulysses data rather than the Monte Carlo simu- 
lations. The resulting CR-induced flow modifications 
are small enough that we do not expect significant 
differences in the particle distributions from a com- 
parable test-particle simulation. Similarly, we have 
included the effects of Alfven wave transport, since it 
would presumably be present in the real interplane- 
tary shocks. 

BOEF95 have studied two shocks. For the first 
shock, observed on April 7, 1991, (hereafter BOEF95- 
1) the following properties are assumed: shock veloc- 
ity, Vg = 153 km s~^, sonic Mach number, Mg = 
6.9; Alfvcnic Mach number, Ma = 3.1; upstream 
field strength, B = 30/LtG; upstream particle den- 
sity, n\ = 1.756cm~^; upstream ion temperature, 
Ti = 3.57 X lO-^K, and magnetic obliquity, 9i = 77°. 
The second shock, which was observed on April 28, 
1991 (hereafter BOEF95-2) was a bit weaker than 
the first shock. This shock is initiated with these 
conditions: shock velocity, Vg = 165 km s~^; sonic 
Mach number, Mg = 3.9; Alfvenic Mach number. 
Ma = 2.2; upstream field strength, B = 20/iG; up- 
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stream particle density, ni = 0.338cm ^, upstream 
ion temperature, Ti — 1.3x lO^K, and magnetic obliq- 
uity, 01 = 75°. 

The grid spacing in all runs for these two shock 
models is set so that Aa; = l/2rg{pth) for the BOEF95- 
1 shock, and Ax = l/^rg{pth) for BOEF95-2, inde- 
pendent of the value of N . Since the shock spreads 
over 3-4 cells in these relatively weak shocks, the effec- 
tive numerical shock thickness, 5x ^ {1 — 2)rg{pth)- 
These values of Ax are necessary to produce parti- 
cle fluxes matching the observations. For Ax twice 
these values, for example, particle fluxes are too low 
to match the observations with any reasonable choices 
of Ci . We note below that the best fits to the Ulysses 
data correspond to iV = 4 for BOEF95-1 and TV = 9 
for BOEF95-2. These shocks are quasi-perpendicular, 
with Oi « 75°and 62 ~ 85°, so that as a particle 
streams a distance A along a field line, it moves along 
the shock normal a distance A cos 9 ^ N cos Org which 
is ^ (1 — 2)rg. Thus, it makes sense that the required 
numerical shock thickness, (5x ~ (1 — 2)rg{pth)- 

The value of p2 was fixed at 3pth,i, where pth,i is 
the momentum at which the Maxwellian distribution 
peaks in the downstream region for the initial shock. 
The CR pressure is small compared to the gas pres- 
sure, and consequently not very different from test- 
particle conditions. Thus, the choice of p2 substan- 
tially affects neither the flow dynamics nor the parti- 
cle spectrum. 

In order to transform the particle distribution func- 
tion (which is isotropic to lowest order in the local 
fluid frame) to the particle count rate in the space 
craft frame, we need to know the velocity of the down- 
stream flow relative to the spacecraft. That velocity is 
difficult to compute accurately from the information 
available, so we chose it to match the particle veloc- 
ity, Vpeak — 500km s~-^at the peak of the MaxweUian 
distribution. 

Fig. 4 shows the computed and measured omni- 
directional particle flux in the spacecraft frame di- 
vided by the particle momentum cubed, p^, and also 
the computed particle distribution function for the 
BOEF95-1 shock. The filled dots are the Ulysses 
data taken from Fig. 1 of BOEF95. Our results 
are shown at t = 6 minutes. For the velocity range 
500 < Vp < 2000 km s^^, the simulated particle 
flux has reached nearly steady values from an ini- 
tial Maxwellian form after 5 minutes. Three values 
oi N = X/vg = 4, 20 and 40 were tried while keeping 
ci = 1.6. For the fourth run, = 4 and ci = 2.0 



were chosen. The grid spacing is the same and so the 
shock thickness is about the same for all four cases. 
All except the = 40 run produce acceptable fits 
to the Ulysses data, although the iV = 4 is some- 
what the best. That value of N was also preferred by 
BOEF95 from their Monte Carlo simulations. The 
similar comparison for BOEF95-2 shock is presented 
in Fig. 5. The same quantities are plotted as in Fig. 
4. Now our results correspond to a shock evolution 
time, i = 10 minutes. The value of ci for the best fit 
is again 1.6. Three values of N are compared; namely 
TV = 9, 20 and 40. For a fourth run iV = 9 and 
ci = 2.0 were used. As in BOEF95, the simulated par- 
ticle fiuxes seem to agree best with observations when 
small values of N are used. Although these calcula- 
tions include full dynamic effects of CRs, the modifi- 
cation to the fiow structure is insignificant as shown 
in Fig. 6. But the slight reduction in postshock pres- 
sure and temperature in dynamic calculations means 
somewhat smaller injection rate compared to the test- 
particle simulations. For test-particle simulations, the 
particle flux shown in Figs. 4-5 could be about 50 % 
larger than that of dynamic runs for the velocity near 
1000 km s^^, for example. 

While, in both examples, the comparisons of each 
case with the BOEF95 fluxes are fairly similar, we can 
see in the p^/(p) plots that smaller N leads to higher 
momentum particles at a given time. That is simply 
due to the fact that smaller N leads to smaller k (see 
equation 3-2), and consequently a smaller acceleration 
time, since the individual particle acceleration time, 
tacc oc K (see, e.g., Lagage & Cesarsky 198S). The 



particle flux near and above pi, however, increases 
with the values of N for three runs with ci = 1.6. 
The particles have larger mean free paths for larger 
N and so have higher probability to cross the shock, 
since the shock thickness is about the same for all 
runs. This leads to a higher particle flux in the in- 
jection pool and so a higher injection rate. However, 
the sensitivity to N is rather weak in our model, com- 
pared to that to ci, for a given shock thickness and 
for a given value of ci. The ci — 2.0 cases produce 
fewer CRs, but accelerate them to the same momenta 
as the same N and ci = 1.6 cases. That is because 
ci = 2.0 places the transition from thermal to non- 
thermal particles farther into the Maxwellian tail of 
the postshock distribution, and thus, reduces the pop- 
ulation of the injection pool. This shows that the 
injection rate is mostly controlled by the choice of ci 
for a given shock thickness. In our numerical injection 
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model we do not have a self-consistent way to deter- 
mine the best value of ci for a given value of N, while 
in Monte Carlo simulations the injection is treated 
self- consistently. On the other hand, the fact that the 
numerical shock thickness must be relatively thin to 
produce consistent fluxes {i.e., Sx ~ {l — 2)rg{pth), so 
that N cos 9 1 for the best fits with N) could imply 
that the observed particle flux cannot be explained if 
the scattering turbulence is weak {e.g., N ^ 1). This 
is consistent with the conclusions of BOEF95. 

Again these comparison calculations have shown 
that the diffusion-convection formalism with our new 
injection scheme and a reasonable set of scattering 
and injection parameters can reproduce the particle 
injection and acceleration processes in real oblique 
MHD shocks. The detailed dependence of our cal- 
culations upon the model parameters such as ci and 
grid spacing should not be overemphasized, since our 
model is not intended to represent the detailed mi- 
crophysics of the injection and shock formation pro- 
cesses, but rather only to try to capture the outcomes 
reasonably well. 

4.3. Tvi^o-Fluid Comparisons 

Beginning from the above successes, it is useful 
to provide direct comparisons between the diffusion- 
convection simulations and the simpler two-fluid ver- 
sions of them. Two-fluid methods have been es- 
pecially useful in complex time dependent applica- 
tions, such as the evolution of supernova remnants 
{e.g., Dorfi 1991; Jones & Kang 1992). They are cur- 



rently the only practical method of calculating multi- 
dimensional CR-modified flows {e.g., ). As mentioned 
in the introduction, there has been some controversy 
in the past about the conditions under which two- 
fluid methods can provide reliable dynamical solu- 
tions for diffusive shock structures. Paper 1 addressed 
some of these issues in the context of parallel shocks, 
and identifies some of the background literature. We 
demonstrated there the basic agreement between two- 
fluid and diffusion-convection methods. Arguments 
are sometimes expressed that momentum-dependent, 
cross-field diffusion in oblique shocks could invalidate 
the fluid-like CR behaviors implicit in the two-fluid 
formalism. To the best of our knowledge the only pre- 
vious comparisons of the methods for oblique MHD 
shocks were by Frank, Jones and Ryu (1995). They 
considered only a case with a momentum independent 
diffusion coefficient and one with weak momentum de- 
pendence, K cx p^/^. Thus, we provide here a similar 



comparison as in Paper I, but now for oblique MHD 
shocks. For this we choose two representative shocks; 
namely, EBJ95-D and BOEF95-1 described in §4.1 
and §4.2, respectively. 

In the two- fluid version of the diffusive acceleration 
model the energy moment of the diffusion-convection 
equation ( [3-lD is integrated from p2 to p3 to produce 
the conservation equation for CR energy; namely, 



dEc 
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where E^, is the CR energy density. No new approxi- 
mations are introduced in deriving equation [4-2| from 
equation 3-1. It does contain three closure param- 
eters, 7c, Stf and (k) that are really properties of 
the solution, but in practice must be estimated a pri- 
ori. For these particular simulations it is sufficient to 
set the CR adiabatic index, 7c = |, since the parti- 
cle populations are entirely nonrelativistic. The in- 
jection energy rate, Stf represents energy exchange 
with the thermal plasma (see also equation [2.5] in 
Paper I). In the diffusion-convection simulations, we 
calculate numerically an analogous injection energy 
rate, S, from the particle flux crossing the low mo- 
mentum boundary of the CR population at p2, and 
subtract it from the thermal energy. Then the spa- 
tially integrated injection rate for the two- fluid model, 
/ = / Sdx, can be parameterized by a dimensionless 
two-fluid "injection parameter" , 77, through the rela- 
tion / = (l/2)u|piui 7] where V2 — P2C (see equation 
[2.4] in Paper I). Thus we calculate r] rather than / 
itself as a function of time for each shock modeled 
using the results of the kinetic equation calculations. 
In practice the value of 77 is fairly constant over time 
in the cases we have considered, so very comparable 
two-fluid solutions would be found by assuming a sin- 
gle value in each test; namely 77 « 0.006 for EBJ95-D 
and r] « 0.002 for BOEF95-1. 

To model the mean diffusion coefficient, (k) (see 
equation [2.13] in Paper I), we used the form for nij)) 
discussed in §3 and a simple power-law model for the 
CR distribution function; namely, J{p) ocp^'^, where 
q = 3r/(r — 1) is the standard test-particle index 
expected for diffusive acceleration and r is the ini- 
tial compression ratio for the shock (see Figs. 2, 6). 
We supposed that the CR distribution extended be- 
tween p2, as defined for the full diffusion-convection 
simulation and ps found from the usual relation be- 
tween particle energy and mean acceleration time 
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( Lagage fc Cesarsky 1983 ; EBJ95). In the present 
context that leads to = p2(l + t/TY^"^, where 
T = |k(p2)/(ui Aw), and Au = U2 — ui, for the ini- 
tial shock. In practice we obtained somewhat better 
matches with the diffusion-convection runs by replac- 
ing the factor | by the factor 2; that is, the distribu- 
tion begins to cut off a little below p3. Our results 
presented here use that latter value. 

Two-fluid models are intended only for dynami- 
cal studies, so the appropriate tests are comparisons 
of shock structures computed by the two-fluid model 
and structures computed by diffusion-convection meth- 
ods (or actual shocks, if available). The two-fluid 
shock structure evolution for EBJ95-D is shown by 
the dotted lines in Fig. 2. The agreement with 
the diffusion-convection simulation is good, reinforc- 
ing our earlier conclusions about the basic validity 
of the two-fluid model. At the end of this simulation 
{t = 1.2x10^* sec) Pc is definitely producing important 
modifications to the flow structure. It is still increas- 
ing, so that a more major dynamical influence could 
be expected at later times. In fact, the time asymp- 
totic two-fluid solution for this shock should be com- 
pletely CR dominated, as can be demonstrated from 
comparable shock solutions in the Figure 3 presented 
by Kang & Jones (1990), or Figure 3 of Frank, Jones 
& Ryu (1995), for example. It is simple to demon- 
strate that the time-asymptotic two-fluid jump con- 
ditions (hence, Pc downstream), arc independent of 
the value or spatial structure of the diffusion coeffi- 
cient. For this shock, however, the time required to 
approach that solution from the initial conditions we 
used would be extremely long; so long, in fact, that 
the practical significance of the time-asymptotic so- 
lution is doubtful. 

By contrast, the shock structure at intermediate 
times is influenced sensitively by the early time evo- 
lution of (k). That grows quickly and asymptotes to 
(k) oc (t/r)'^. The rate of dynamical shock evolu- 
tion generally scales inversely as the "diffusion time" , 
td = That means at late times the shock 

structure evolves very slowly. At intermediate times, 
Pc is largely controlled by the very early history of the 
shock; particularly {k) and "q. So, the match we see in 
Fig. 2 between the diffusion-convection and two-ffuid 
simulations is sensitive to the values of r and k{p2). 
That the appropriate value of r is reasonably well pre- 
dicted by the simple test particle theory for the diffu- 
sion coefficient is an encouraging outcome. The mi- 
nor differences between the diffusion-convection and 



two-fluid runs come, in fact, from small differences in 
the upstream spatial variations of (k) modeled in the 
two-fluid calculations and as computed directly from 
the diffusion-convection simulation. In the diffusion- 
convection simulation, the particle distribution tends 
to "harden" signiflcantly upstream of the shock (see, 
e.g., Kang & Jones 1991; Paper I), so the actual (k) 
increases away from the shock, reducing the value of 
Pc as a result. 

The BOEF95-1 shock has a dynamically very sig- 
niflcant magnetic field, so it presents an important 
comparison case for two-fluid models in the MHD 
shock context. Both the two-fluid and the diffusion- 
convection shock structures for BOEF95-1 are shown 
in Fig. 6. The two-ffuid parameters were deter- 
mined in exactly the same way as for EBJ95-D. Again 
the agreement between the two models is very good. 
Recall that we already made a comparison between 
the particle velocity distributions from the diffusion- 
convection solution and the Ulysses observations, but 
that we have no detailed information about the phys- 
ical, interplanetary shock structure. In this case the 
value of Pc at the shock is less than 10% of the gas 
pressure, Pg , by the end of the simulations (< = 6 min- 
utes), so there is only minor modification of the shock 
structure by the CRs. In BOEF95-1, as in EBJ95-D 
the structure at intermediate times is primarily de- 
termined by the evolution of (k) at the early times in 
the simulation. Again, because (k) increases rapidly 
with time, the td that determines the rate of shock 
evolution becomes very long. A direct consequence 
of this is that, despite the apparently steady shock 
by the end of the simulation (and also the diffusion- 
convection simulation), we are not seeing the solution 
that would be determined from the steady-state two- 
fiuid model. 

To illustrate the point, we can take advantage 
of the argument made earlier that the steady-state 
jump conditions can be found using any convenient 
diffusion coefficient. For this we have recomputed 
BOEF95-1 using a constant {n) = 0.4, and allowed 
it to evolve to a steady state. In this case we use the 
criterion for a time-asymptotic solution that not only 
does the peak value of Pc become steady, but also that 
Pc be uniform directly behind the shock. The evolu- 
tion of this shock {BOEF95-1C) and its final struc- 
ture arc shown in Fig. 7. Early on, the evolution of 
BOEF95-1C and BOEF95-1 are similar, because (k) 
are comparable. However, as the BOEF95-1 diffu- 
sion coefficient increases dynamical evolution "stalls" , 
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while the constant diffusion coefficient of BOEF95- 
IC allows it to continue directly towards the formal 
steady state solution. The final acceleration efficiency 
of this shock as measured by Pc is more than an order 
of magnitude greater than for the BOEF95-1 simula- 
tions presented earlier or the physical shock. Thus, 
it would clearly be inappropriate to apply the steady 
state two-fluid model to estimate the efficiency ex- 
pected in this shock for even moderately long times. 

The discrepancy is not an indication of basic flaws 
in the two-fluid model, but rather that time-asymptotic 
solutions are not very relevant in this situation. The 
key question becomes how to choose, in the two-fluid 
model, a meaningful set of closure parameters for 
a time dependent calculation. In both of the tests 
conducted here it seems adequate to apply simple 
test-particle models to those parameters, because the 
shock structure is not sufficiently modified on short 
times. Over longer times that convenience may be- 
come dubious, but, since the shock properties at mod- 
erately late times are largely set by the conditions 
early in the shock evolution, this breakdown may not 
be crucial in many instances. It is important to know, 
then, if the shock structure should evolve quickly on 
time scales of interest. That is something wc can hope 
to estimate reasonably well using the standard test- 
particle approach. These issues are more important if 
we expect a steep momentum dependence to the diffu- 
sion coefficient, as in the cases studied here, so a more 
basic understanding of the evolution of the resonant 
Alfven wave field becomes important, as well. 

5. Conclusion 

In order to study the particle acceleration in oblique 
magnetohydrodynamic (MHD) shocks, we have im- 
plemented the existing diffusion-convection methods 
of Kang & Jones 1991 into a full MHD code, and 
adopted a "thermal leakage" type injection model in- 
troduced by Kang & Jones 1995 (Paper I). In our 
injection model, the distribution of the suprathermal 
particles which cannot be treated properly with the 
diffusion-convection method was assumed to match 
smoothly onto the Maxwellian distribution of the 
gas particles. The matching condition is controlled 
by a free parameter ci, which in turn determines 
the particle injection rate into the CR population. 
Firstly, we have calculated the MHD shocks for var- 
ious field obliquities considered by Ellison, Baring & 
Jones (1995), in order to study the dependence of the 



injection efficiency on some shock properties via test- 
particle Monte Carlo simulations. By adjusting the 
free parameter ci of our injection model over a modest 
range wc were able to demonstrate that our numer- 
ical technique can, in fact, produce particle spectra 
comparable to theirs. Secondly, we have reproduced 
the proton flux distributions at oblique interplane- 
tary shocks observed in situ by the Ulysses spacecraft. 
These shocks have also been simulated previously by 
Baring et al., 1995, via test-particle Monte Carlo tech- 
nique. We adopted the shock parameters chosen by 
them to match direct observations. To obtain good 
agreement with the observations, the numerical shock 
thickness for these quasi-perpendicular shocks has to 
be about (1-2) times the thermal gyro-radius. This is 
consistent with the conclusion of Baring et al., that 
strong scattering turbulence was present in these in- 
terplanetary shocks. 

The Monte Carlo technique treats both thermal 
and cosmic ray particles by the same scattering law, 
so the injection process comes about naturally via 
the acceleration of thermal particles to higher en- 
ergies. In contrast, injection cannot be determined 
self-consistently through diffusion-convection models 
for cosmic-ray transport, since particles at momenta 
where injection occurs do not form an isotropic distri- 
bution and the diffusion approximation is not valid. 
Our model works around this by introducing a free 
parameter that establishes the momentum at which 
the suprathermal distribution must match onto the 
thermal distribution behind the shock. One might 
be concerned that the particle distribution at inter- 
mediate energies (between the thermal peak energy 
and energies much greater than thermal peak energy) 
would be dependent upon the details of the injection 
process. Our simulations, however, seem to indicate 
that the dependence is not sensitive enough to make 
a clear distinction between the particle spectra simu- 
lated with our diffusive injection model and existing 
observations, or the particle spectra simulated with 
Monte Carlo techniques. This leads us to the tenta- 
tive, but encouraging conclusion that a simple, macro- 
physical model like ours can offer a practical compro- 
mise between ad hoc injection and injection models 
built to include microphysical details. 

Comparison tests presented here and in Paper I 
have shown that, in our model, the injection pro- 
cess and its rate are mostly determined by the nu- 
merical shock thickness in terms of the thermal gyro- 
radius and the momentum, pi = c\ pth, where the 
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suprathermal particle' distribTition matches onto the 
Maxwellian distribution. Presumably both of these 
are affected in a detailed model by the strength of 
Alfven turbulence {e.g., Malkov & Volk 1995). Under 
the assumption that the shock thickness is about the 
mean scattering length of the thermal peak momen- 
tum, pth, for quasi-parallel shocks and the thermal 
gyro radius, rg{pth) for quasi-perpendicular shocks, 
the appropriate value of c\ lies between 1.5-2. A 
smaller ci leads to larger injection rates, because it 
allows a larger "injection pool" of diffusive particles 
to form. EBJ95 found in Monte Carlo simulations 
that the rate of particle injection decreased strongly 
with increasing obliquity for strong shocks, unless the 
cross-field diffusion was strong. In our model that is 
effected by increasing the parameter ci with obliq- 
uity or decreasing it towards smaller Mach number or 
stronger Alfven turbulence. More quantitative pre- 
diction requires much further work, however. We note 
also that the injection efficiencies given in EBJ95, 
which were based on test-particle simulations, are 
overestimated compared to dynamical shocks, espe- 
cially for strong, quasi-parallel shocks, since the CR 
energy density predicted by those simulations is sig- 
nificant enough that the CRs would modify the shock 
dynamics. In particular, the thermal particle dis- 
tribution would be colder and the subshock velocity 
jump would be smaller in self-consistent dynamic cal- 
culations. These differences between test-particle and 
dynamical shocks were confirmed by the recent Monte 
Carlo simulations of EBJ96. 

In Paper I, we showed that diffusive acceleration 
numerical methods applied to parallel shocks pro- 
duce similar shock structures and particle distribu- 
tions compared to Monte Carlo and hybrid plasma 
methods, and that they matched direct observations 
at the earth's bow shock. The comparisons reported 
here for oblique shocks strengthen the important con- 
clusions of Paper I that the essential physics of the 
particle injection and acceleration can be captured 
by each of these diverse computational methods, and 
that they are all practical and complementary tools 
for understanding the physics of diffusive shock ac- 
celeration. This also implies that our numerical ap- 
proach provides a way to incorporate naturally the in- 
jection process into the existing diffusion-convection 
technique. The advantages of this formalism distin- 
guishing it from Monte Carlo or plasma methods are 
that it is time-dependent, in addition to being a fully 
dynamically self-consistent MHD diffusion-convection 



tecimique, so that it can be used for evolving and 
structurally complex flows. In upcoming papers, we 
will use it to study the acceleration efficiency and the 
nonlinear back reaction of CR pressure on the shock 
dynamics in various astrophysical shock waves. 

We also simulated a pair of oblique two- fluid shocks. 
Each was constructed exactly as for one of the diffusion- 
convection simulations reported, with the required 
closure parameters determined from simple test-particle 
considerations. The dynamical properties of the two- 
fluid shocks are quite consistent with the diffusion- 
convection solutions. These simulations demonstrate 
in the oblique shock context the basic validity of 
the two-fluid method. We emphasize, however, that 
steady state two-fluid solutions may not be applicable, 
even when the shock structures appear to be steady. 
If the cosmic-ray diffusion coefficient has a strong 
momentum dependence, the rate of shock evolution 
can become very slow, so that while a shock may ap- 
pear dynamically steady, in practical terms the time- 
asymptotic solution is not likely to be reached for a 
long time. Then the dynamical conditions creating 
the shock may very likely have changed, requiring the 
shock to readjust once more. 

We are grateful to Matthew Baring for illuminat- 
ing discussions about oblique CR shocks and Monte 
Carlo methods and for useful suggestions that helped 
us improve the manuscript. This work was sup- 
ported in part at the University of Minnesota by the 
NSF through grant AST-9318959, by NASA through 
grants NAGW-2548 and NAG5-50505 and by the Uni- 
versity of Minnesota Supercomputer Institute. HK is 
supported in part at Pusan National University by 
the Korean Research Foundation through the Brain 
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Fig. 1. — Top panel shows the distribution function g = f{xs,p)p'^ versus the particle momentum p in units of 
mc for particles at the shock for EBJ95 test runs. The lines are labeled by the value of the upstream obliquity 
angle 9i. The results are shown at (t/lO^s) = 12,8,6, and 4 for 01 = 0°, 20°, 30°and 35°, respectively. See text 
for the shock parameters. Bottom panel shows the integral density distributions calculated from the momentum 
distribution function shown in the top panel. The open and filled circles and squares are representative points of 
EBJ95 Monte Carlo simulation results of the same shock conditions (from their Fig. 5). 
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Fig. 2. — The shock flow structure of the model EBJ95-D shock at (t/lO^s) = 0.4, 0.8, 1.2 in our fuUy dynamical 
simulations (solid line). The initial condition is specified by the pure MHD shock jump. The two-fluid solution 
for the same shock as discussed in §4.3 is shown by the dotted line. This is a Mach 100 shock, with a very weak 
magnetic field at an upstream obliquity, 9\ = 30°. 
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Fig. 3. — Top panel shows the distribution function g = f{xs,p)p'^ versus the particle momentum p at the shock 
for EBJ95 test runs at (t/lO^s) = 1.2. The obHquity, 9i = 30° and the injection parameter, ci = 2.0 The solid line 
is the spectrum from the test-particle run, while the dotted of line is for fully dynamic runs. Bottom panel shows 
the integral density distributions calculated from the distribution functions shown in the top panel. 
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Fig. 4. Simulated omni-directional particle flux in the Ulysses spacecraft frame divided by the particle momentum 
cubed, (upper), and the particle distribution function f{p) (lower) for BOEF95-1 shock at t = 6 minutes. The 
solid line is for = 4, dashed line for N = 20, and long-dashed line for N = 40. The value of Ci = 1.6 for these 
three runs. The dotted line is for A'' = 4 and ci = 2.0 
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Fig. 5. Same as Fig. 4, except for BOEF95-2 shock at t = 10 minutes. The sohd hue is for N = 9, dashed Hne 
for N = 20, and long-dashed hne for N = 40. The value of Ci = 1.6 for these three runs. The dotted line is for 
A'' = 9 and ci = 2.0 
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Fig. 6. — The shock structure computed for the BOEF95-1 shock at t = 2, 4, 6 minutes. The diffusion-convection 
solution is shown by the sohd hues and the two-fluid solution by the dotted lines. The two-fluid solution uses a 
mean diffusion coefficient that evolves in time according to a simple test particle model for the distribution function. 
The shock initial conditions are indicated by the discontinuous curves. Details are given in the text. 
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Fig. 7. — Evolution of the shock structure for the shock BOEF95-1C. This is a two-fluid model shock and differs 
from the two-fluid shock shown in Fig. 6 only in the use here of a constant diffusion coefficient, (k) = 0.4. Times 
represented are f = 0, 6, 24, 72. For the last time the shock has approached its time-asymptotic limit. 
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